!> \brief DE-RSM interface.
!!
!! \details This module contains subroutines related to the hybridization
!! of the Differential Evolution algorithm with the Response Surface Methodology.
module hybrid

   use rsm, only: get_optimum_solution

   implicit none

contains


   !> \brief Checks if the population size is enough to apply the RSM
   subroutine checks_rsm_condition( kh, nu, np, nf)
      implicit none
      integer, intent(in)  :: kh    !< Kind of hybridization ( 0 = none, 1 = RSM-min, 2=RSM-max)
      integer, intent(in)  :: nu    !< Number of unknowns
      integer, intent(in)  :: np    !< Size of the population
      integer, intent(out) :: nf    !< Number of individuals for fitting

      if ( np < (nu+1)*(nu+2) / 2 ) then

         write(*,*) "ERROR: population size insufficient to apply RSM."

         stop

      end if

      if ( kh == 1 ) then

         nf = (nu+1)*(nu+2) / 2

      else if ( kh == 2 ) then

         nf = np

      else

         write(*,*) "ERROR: unknown kind of hybridization."

         stop

      end if

   end subroutine checks_rsm_condition



   !> \brief Returns the best estimate given by the Response Surface Methodology
   subroutine get_rsm_optimum( nf, nu, np, pop, fit, x, es)
      implicit none
      integer, intent(in) :: nf !< Number of selected individuals for fitting
      integer, intent(in) :: nu !< Dimension of the problem
      integer, intent(in) :: np !< Size of the population
      real(8), dimension(np,nu), intent(in)  :: pop !< Population
      real(8), dimension(np),    intent(in)  :: fit !< Fitness of the population
      real(8), dimension(nu),    intent(out) ::   x !< Best estimated individual
      integer,                   intent(out) ::  es !< Exit status: 0 = success, 1 = failure


      ! Inner variables

      integer ::   i, j ! Dummy index
      integer ::    nf1 ! Variable nf
      integer ::  ibest ! Index of the best individual
      integer :: iworst ! Index of the worst individual

      real(8) ::   rmin ! Auxiliary variable
      real(8) ::   rmax ! Auxiliary variable

      real(8), dimension(np) ::  dist ! Distance to the best individual
      integer, dimension(np) ::    od ! Ordering of vector dist

      real(8), dimension(np, nu) :: dm ! Design matrix
      real(8), dimension(np)     :: fm ! f

      ! Searching for the best and worst individual of the population

      ibest = 1

      iworst = 1

      rmin = fit(ibest)

      rmax = fit(iworst)

      do i = 1, np

         ! Minimum
         if ( fit(i) < rmin ) then

            rmin = fit(i)

            iworst = i

         end if

         ! Maximum
         if ( rmax < fit(i) ) then

            rmax = fit(i)

            ibest = i

         end if

      end do

      ! Calculating the distance of each individual to the best one

      do j = 1, np

         ! Temporary use of x (to reduce memory and CPU time)

         x = pop(j,:) - pop(ibest,:)

         dist(j) = sqrt( dot_product( x, x ) )

      end do

      ! Ordering 'dist'

      do j = 1, np

         od(j) = j

      end do

      call qsort2(dist,od)


      ! Ordering the individuals closer to the best one

      do j = 1, np

         dm(j,:) = pop(od(j),:)

         fm(j) = fit(od(j))

      end do

      ! Searching for the best solution
      ! ( The following code tries to apply RSM while the number of
      !   fitting functions is less than or equal to the size of the
      !   population or a success is returned by the RSM subroutine )

      es = 1

      nf1 = nf

      do while ( es /= 0 .and. nf1 <= np )

         call get_optimum_solution(nf1, nu, dm(1:nf1,:), fm(1:nf1), x, es)

         nf1 = nf1 + 1

      end do

   end subroutine get_rsm_optimum


   !> \brief Sorts a vector x based on the Quicksort Algorithm.
   !! It also returns the ordering ox of the original vector.
   recursive subroutine qsort2(x,ox)
      implicit none
      real(8), dimension(:), intent(inout) ::  x !< Unordered vector
      integer, dimension(:), intent(inout) :: ox !< Ordering of x

      ! Inner variables

      integer :: ip ! Split position

      if ( size(x) > 1 ) then

         call partition2( x, ip, ox )

         call qsort2( x(:ip-1), ox(:ip-1) )

         call qsort2( x(ip:), ox(ip:) )

      end if

   end subroutine qsort2


   !> \brief Separates the number lower and greater than the pivot
   subroutine partition2(x, ip, ox)
      implicit none
      real(8), dimension(:), intent(inout) :: x  !< Unordered vector
      integer,                 intent(out) :: ip !< Split position
      integer, dimension(:), intent(inout) :: ox !< Ordering of x

      ! Inner variables

      integer :: left, right, iaux

      real(8) :: pivot, raux

      ! Initializing sweepers

      left = 0

      right = size(x) + 1


      ! Pivot

      pivot = ( x(1) + x(size(x)) ) / 2


      ! Sweeping unordered vector: right -> left and left -> right

      do while ( left < right )


         ! Sweeping unordered vector: right -> left

         right = right - 1

         do while ( pivot < x(right) )

            right = right - 1

         end do


         ! Sweeping unordered vector: left -> right

         left = left + 1

         do while ( x(left) < pivot )

            left = left + 1

         end do


         if ( left < right ) then

            raux = x(left)

            x(left) = x(right)

            x(right) = raux

            iaux = ox(left)

            ox(left) = ox(right)

            ox(right) = iaux

         end if

      end do

      if ( left == right ) then

         ip = left + 1

      else

         ip = left

      end if

   end subroutine partition2


   !> \brief Generates a mean individual
   subroutine get_mean_individual(nu, np, pop, x)
      implicit none
      integer, intent(in) :: nu !< Dimension of the problem
      integer, intent(in) :: np !< Size of the population
      real(8), dimension(np,nu), intent(in)  :: pop !< Population
      real(8), dimension(nu),    intent(out) ::   x !< Mean individual

      ! Inner variables

      integer :: i      ! Dummy index


      ! Calculating the mean individual

      x = 0.d0

      do i = 1, np

         x = x + pop(i,:)

      end do

      x = x / dble(np)

   end subroutine get_mean_individual


end module hybrid
